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ABSTRACT 

The intermediate phases of planet formation are not directly observable due to lack of emission 
from planetesimals. Planet formation is, however, a dynamically active process resulting in collisions 
between the evolving planetesimals and the production of dust. Thus, indirect observation of planet 
formation may indeed be possible in the near future. In this paper we present synthetic observations 
based on numerical N-body simulations of the intermediate phase of planet formation including a 
state-of-the-art collision model, EDACM, which allows multiple collision outcomes, such as, accretion, 
erosion, and bouncing events. We show that the formation of planetary embryos may be indirectly 
observable by a fully functioning ALMA telescope if the surface area involved in planetesimal evolution 
is sufficiently large and/or the amount of dust produced in the collisions is sufficiently high in mass. 
Subject headings: planets and satellites: formation, methods: numerical 


1. INTRODUCTION 


Planet formation produces a broad range of plane¬ 
tary outcomes ranging fro m hot Jupiters around main- 
sequence so lar-type stars (Mayor & Queloz 1995 Marcy 


et al.|1997 ) to lunar- mass planets around evolved degen¬ 


erate" neutron stars (|Wolszczan & Frail| 1992). Planets 
appear to be a common by-product of star formation 
with at least one i n every four adult sta rs found to have at 
least one planet ( Petigura et al.| [2013). Recent advances 
in observational techniques and hardware have led to a 
huge inc rease in the number and diversity of extrasolar 
planets (Lissau er et al.|[2014 > and unprecedented resolu¬ 
tion and detail in young protoplanetary disks. However, 
these observations represent snapshots or still frames of 
the beginning and end of a dynamic evolving story of 
planet formation giving us just the briefest look behind 
the scenes. 

The intermediate stages of planet formation are diffi¬ 
cult to observe directly for a number of reasons: first, 
planetesimals, the building blocks of planets, have a 
small surface area to mass ratio thus they are difficult 
to observe in reflected light from the star, second, plan¬ 
etesimals are relatively low-mass objects, thus, they do 
not produce significant thermal energy themselves nor do 
they produce any significant gravitational effect on the 
central star; third, young stars are notoriously noisy at a 
number of wavelengths making searches for faint objects 
near the central star more difficult. In the core accre¬ 
tion model of planet formation planetesimals evolve into 
planets and embryos via collisions with one another. All 
collisions whether accretion dominated or erosive result 
in the production of some dust and dust has a high sur¬ 
face area to mass ratio making it highly detectible at 
infra-red and sub-mm wavelengths. Thus, although di¬ 
rect observation of planetesimal evolution and planet for¬ 
mation may not be possible in the near future it may be 
more prudent to search for indirect signatures by looking 
for collisionally generated dust. 

In this work we present numerical simulations of terres¬ 
trial planet formation including a state-of-the-art plan¬ 
etesimal collision model, which can describe a range of 


collision outcomes from merging to catastrophic disrup¬ 
tion. Using the collisions from the evolving planetesimal 
population we construct synthetic dust images and find 
that terrestrial planet formation may just be observable 
with a fully operational ALMA telescope. 

2. NUMERICAL METHOD 


All simulations were carried out using the parallelised 


IV-body gravity code PKDGRAV (Stadel 2001 

Richard- 

son et al. 

2000 Leinhardt & Richardson 2005 

). PKD- 


GRAV is a second order leap-frog integrator with a hi¬ 
erarchical tree to efficiently calculate gravitational inter¬ 
actions between large numbers of particles. 


2.1. Planetesimal Disk Model 

All simulations presented here are of modest resolu¬ 
tion (N = 10 4 ) in order to compare directly with previ¬ 
ous work (|Kokubo & Ida|[2002| |Leinhardt & Richardson 
2005|). The simulations begin with equal-sized particles 
arranged in an annulus from 0.5 to 1.5 AU around a 
single central stellar potential of one solar mass. The 
particles are distributed with a power-law surface den¬ 
sity distribution E = Ei(a/AU) _ “, where a = 3/2. The 
total mass in solid material is 2.8M® similar to the mass 
derived from the minimum mass solar nebula. 

Planetesimal particles were radially expanded by a fac¬ 
tor / in order to accelerate the planetesimal evolution. 
The inflation of the planetesimal particl es reduces the 
evolu tion timescale by approximately / 2 (Kokubo & Ida 
2002). The specific value of the expansion factor (/ = 6) 
was chosen to be consistent with previous work. I t has 
been shown analytically in Kokubo & Ida (1996) that 
/ ~ 6 should not significantly impact the growth phase 
of planetesimals as long as gravitational scattering is rel¬ 
atively unimportant and the planetesimals are evolving 
in a simple system with a single central potential and 
no giant planet perturbers. For simplicity nebular gas 
is neglected in this work. Planetesimals will be negli¬ 
gibly effected by this simplification (they are large and 
would suffer little aerodynamic drag though collisionally 
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generated debris would be more substantially affected) Q 

2.2. Planetesimal Collision Model 

In this paper we compare terrestrial planet forma¬ 
tion simulations us ing t hree c oll ision models: p erfect 
merging, RUBBLE ([Leinhardt & Ri chardson 2005) , and 
our n ew collision model EDACM (|Leinhardt &; Stewart 
20121. In both the RUBBLE and EDACiVl collision mod¬ 
els two timesteps are used within the simulation to allow 
particle collisions to be accurately resolved but keep the 
simulation as efficient as possible. The major “orbital” 
timestep is determined by the orbital dynamical time 
(~ 1 yr at 1 AU), however, the “collisional” timestep 
is determined by the planetesimal dynamical time (~ 
1 /yjGp) which is of order hours (p = 0.00925gem -3 for 
/ = 6). Therefore, in all EDACM and RUBBLE simula¬ 
tions all planetesimals were initially placed on the major 
timestep (0.01 yr), when a collision was detected the im- 
pactors were demoted to the minor collisional timestep 
(1.5 x 10~ 4 yr) and the gravity is calculated about 64 
times for the colliding particles and once f or all other 


particles (for detai ls see Leinhardt & Richardson 2005 
B onsor etal.| |2015). 


2.2.1. Perfect Merging 

The simplest collision model used in these simulations 
is perfect merging. In our implementation of this model 
two colliding planetesimals merge on impact assuming a 
perfectly inelastic collision if the impact speed is less than 
or equal to the mutual particle escape speed regardless of 
impact angle or mass ratio. If the impact speed is greater 
than the particle escape speed the impactors bounce in- 
elastically with a normal coefficient of restitution of 0.8. 
This collision model is the most common model used in 


previous IV-body sim ulations (for example, Kokubo & 
Ida||1998[ [2000l [2002] ). 


2.2.2. RUBBLE 

In addition to perfect merging we compare our new 
collision model to RUBBLE the collision model used in 


Leinhardt & Richardson (2005). RUBBLE is nominally 
more accurate than perfect merging and allows for a va¬ 
riety of collision outcomes, however, it does assume that 
collisions between planetesimals are subsonic and cause 
minimal deformation and phase changes. In RUBBLE all 
planetesimals a re assumed to be gravit ational aggregates 
or rubble piles (Richardson et al. 20021. When a collision 
is detected the outcome is determined in the first instance 
by using a look-up table of pre-calculated rubble-pile col¬ 
lisions. If the largest remnant mass is greater than 80% 
of the total colliding mass the collision is deemed “sim¬ 
ple” and the second largest remnant is small, the largest 
remnant from the look-up table is used as the collision 
outcome and any remnant mass is relegated to unresolved 
debris. If the collision is “complex”, the look-up table re¬ 
turns a largest remnant mass that is less than 80% of the 
total colliding mass, the planetesimal particles are sub¬ 
stituted with rubble-piles (gravitational aggregates of in¬ 
destructible billiard balls) and the collision is calculated 

1 We are currently working on simulations with a lower ex¬ 
pansion factor and including nebular gas that would allow us to 
complete more diverse simulations that evolved planetesimals to a 
later state (see Carter et al. in prep). 



Figure 1. Snapshots at t = 0,1, 2, 3 sec for frames A-D showing 
an example of a “perfect” hit-and-run collision integrated using the 
iV-body code PKDGRAV and the EDACM collision model. The 
impact shown is a grazing (b = 0.7), equal-mass (particle radius is 
1 km) collision with an impact speed of 50 m s^ 1 . 


directly within the planet formation evolution simula¬ 
tion. Once the collision has completed any collisional 
remnants that are less massive than the resolution limit 
are relegated to unresolved debris. Any collisional rem¬ 
nants more massive than the resolution limit continue to 
be followed directly by PKDGRAV. 

The unresolved debris is distributed in the planet for¬ 
mation region in ten circular annuli. Any mass that is 
relegated to the unresolved debris is added to the an¬ 
nuli at the instantaneous location of the collision. The 
unresolved debris is assumed to have circular Keplerian 
orbits. Resolved planetesimals can accrete unresolved 
debris in proportion to their geometric cros s-sectio n and 


the eccentricity o f their orbit (for details see Leinhardt & 


Richardson [2005 ) as they move through the unresolved 
debris annum 


2.2.3. EDACM: Collision Outcomes 

Although the RUBBLE collision model is more accu¬ 
rate than perfect merging, allowing the production of 
debris, RUBBLE becomes less accurate at high impact 
speeds because there is no provision for damage or phase 
changes within the planetesimals. This then limits the 
versatility of the RUBBLE model. In addition, resolving 
the more complex collisions is computationally expensive 
meaning RUBBLE simulations are considerably slower to 
run than their perfect merging counter parts. In order to 
rectify some of these draw backs we developed a new em¬ 


pirically derived analy tic collision model, EDACM (Lein- 
hardt & Stewart|2012 ), which is both more versatile and 
com puta tionall y taster. 


Leinhardt & Stewart 
tion ot the EDACM co 


(2012) presents a detailed deriva- 
lision model and IStewart &; Lein-I 
hardt (|2012|) presents a succinct summary. Thus, here we 
will only briefly describe the model and instead highlight 
differences or expansions to the original model. 

In the PKDGRAV implementation of EDACM there 
are five collision outcome regimes: 


1) Perfect merging - two colliders hit each other 
and merge into one body. This occurs when the impact 
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speed, Vi, is less than the effective escape speed, 

VLc = y/2GM' ot /(R P + Rt ), (1) 

where M[ ot = aM p +M t , and M p is scaled by a, the mass 
fraction of the projectile that geometrically overlaps with 
the target: 

_ J 1 : Rt. > b(R p + Rt) + R p 

— Yp{ttR p 1 2 — 7rZ 3 /3): otherwise, ^ ’ 

where b is the impact parameter, and l/(2R p ) is the 
fraction of the projectile diameter overlapping the 
target at impact, R p , M p , R t , and M t are the radii and 
masses of the projectile, and target, respectively. In this 
collision outcome PKDGRAV replaces the two colliders 
with a single particle conserving mass and momentum. 


2) Hit-and-run - a collision in which the target 
and projectile have a grazing encounter, the target 
looses no mass but the direction of travel is altered. The 
projectile may or may not suffer erosion depending on 
the impact energy. 

A hit-and-run event occurs when the impact speed is 
above V/ sc but below the velocity needed for erosion, 


V er0s i 0 n — \j2M to tQ erosion //b 


( 3 ) 


where M to t is the total system mass, p is the reduced 
mass, and Qerosion is the specific energy necessary to cre¬ 
ate a lar gest remnant that is the ma ss of the target (see 
step 4 in Leinhardt & Stewart||2012j) for the parameters 
of the collision (mass ratio, and impact parameter). In 
addition, for the collision outcome to be in the hit-and- 
run regime the impact must be a grazing event where 
the impact parameter is greater than the critical impact 
parameter, 

bcrit = Did (4) 

K v + Kt 


which is defined by Asphaug (20101 as the impact pa¬ 
rameter where the centre of the projectile is tangent to 
the surface of the target. 

The projectile will remain intact if specific impact en¬ 
ergy from the interacting target mass is less than the 
erosion threshold for the projectile. In PKDGRAV this 
type of collision is considered a perfect hit-and-run colli¬ 
sion and is modelled as an inelastic bouncing event with 
a normal coefficient of restitution of 0.8. An example of 
a “bouncing” hit-and-run is shown in Fig. |T] 

If the impact is more energetic and the specific impact 
energy exceeds the erosion threshold for the projectile, it 
will erode. In this case the largest remnant is the origi¬ 
nal target and the second largest remnant is determined 
using the universal law which is normally written as 


M lr /M tot = -0.5 (Qr/Q^ d - 1) + 0.5, 


( 5 ) 


where M; r is the largest remnant, Qr is the specific im¬ 
pact energy, and Qr D is the critical disruption criterion 
for the specific collision being considered, or 

Mi r /Mtot = ^(Qr/Q'/id) 71 , ( 6 ) 

if Qr > 2 Qrr>, meaning the collision is in the super- 
catastrophic regime, where r) = —1.5 (see appendix 



Figure 2. Snapshots at t = 0,1,1.5, 2 sec for frames A-D, re¬ 
spectively, of a partial erosion collision outcome using EDACM 
implemented PKDGRAV as in Fig. [2] In this impact the colliders 
are again equal-mass with a bulk density of 1 g cm -3 but have 
a larger radius (r = 10 km) than the example shown in Fig. m 
a lower impact parameter (6 = 0.2), and a smaller impact speed 
(Pi = 40 m s -1 ). In frame D the largest remnant, which is 8 km in 
radius, is shown in green and the debris field is shown in grey. The 
viewing angle shown is perpendicular to the debris field, meaning 
this debris field is actually a disk. 


Leinhardt & Stewart 2012). The collisional debris gen¬ 
erated is from the erosion or disruption of the projectile 
only, thus, the Mi r is the largest remnant from the 
projectile disruption and is actually the second largest 
remnant in the collision outcome which includes the 
target. Thus, M to t is simply the mass of the projectile. 
The initial position, veloc ity, an d mass of the debris 
field is described in Section 12.2.61 below. 


3) Partial accretion - two colliders hit each other and the 
target gains mass. A partial accretion collision resolved 
with EDACM would look similar to the partial erosion 
collision shown in Fig. [2] but with the largest remnant 
more massive than the target. This type of collision 
event occurs in a non-grazing collision (b < b cr u ) when 
V' esc < Vi < V eros i on . In this case Mr is determined 
using the universal law Eq. [5] and the distribution of 
the debris field is determinedm the same way as it is 
when a pro jectile is disrupted in a hit-and-run event 
(see Section 2.2.6 for details). 


4) Erosion - two colliders hit each other and the 
target loses mass (see Fig. [2]). A partial erosion event 
occurs when V > V eros i on and 0.1 M t < Mi r < M t . 
In this regime M\ r is still determined by Eq. [5] In 
PKDGRAV this means the collision outcome results 
in the target particle losing mass and the possible 
generation of a debris cloud of resolved particles. 


5) Supercatastrophic Disruption - two colliders im¬ 
pact each other at high impact speeds, large in 
comparison to V' esc , resulting in M\ r < 0.1 M t . The small 
mass of the largest remnant means there is effectively 
no separation in the size distribution between Mj r and 
the second largest remnant, M s i r . In impacts of this 
speed Mi r is determined by the super-catastrophic 
scaling law Eq. [6] The only difference in implementation 
within PKDGRAV is the way that the Mr is deter¬ 
mined. The debris field is generated in the same way 
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as it is in partial accretion and erosion collision outcomes. 


2.2.4. EDACM: Resolution Limit 

Since the EDACM collision model can generate par¬ 
ticles as the result of a collision, similar to RUBBLE, 
we impose a resolution limit, which in most cases is 
the initial planetesimal size. Collisional remnants with 
sizes smaller than the resolution limit are treated semi- 


analytically (see end of Section 2.2.2). Any collisional 
remnants that are larger than the resolution limit are 
followed directly by PKDGRAV. 

2.2.5. EDACM: Post-collision Particle Orbits 

The orbits of the particles after collision depend on b. 
In the case of a perfect merging event the one remnant 
conserves momentum and is placed on the centre of mass 
orbit. In all other collision outcomes, in which the largest 
remnant is resolved, is placed at the centre of mass 
with the centre of mass velocity, V com , if the collision 
is head-on (6 = 0). If the impact has a high impact 
parameter (6 > 0.7) Mi r takes on the velocity of the 
target, V t . If the impact was somewhere in-between (0 < 
6 < 0.7) we assume a linear function for the velocity 
vector, 


V, r = 


Vtb + Vr. 


>(0.7-6) 


0.7 


( 7 ) 


as suggested in Leinhardt & Stewart (2012). 


2.2.6. EDACM: Debris Field 


The size distribution of the debris field (all debris less 
massive than the largest remnant) is determined by as¬ 
suming a differential power-law size distri bution with an 
index of -2.85 (|Leinhardt & Stewart|2012|). The remnant 
masses smaller than and including the second largest 
remnant were determined by the power-law and conser¬ 
vation of mass (total mass in the debris tail could not 
exceed M tot - M ir ). 

The velocity distribution of the debris is determined 
by conserving momentum and ensuring that the total 
energy of all collisional remnants did not exceed the ini¬ 
tial energy of the colliders. Energy is allowed to be lost, 
however. For simplicity energy equipartition amongst 
the debris is assumed meaning that the most massive 
remnant in the debris field will have the slowest speed 
and the least massive will ha ve the highest speed. Al¬ 
though pLiemhardOT^tewart (2012) find that the small 
debris should be found at ail speeds the larger debris does 
travel at the slowest speeds and there was not enough res¬ 
olution in their simulations to determine how the velocity 
should be distributed over the remnants more accurately. 

The spatial distribution of the particles in the debris 
field is complex. In some impacts debris forms a disk 
around the largest remnant but in many collisions the de¬ 
bris forms two jet-like features (Fig. [3]). In some instances 
these “jets” are imbalanced and remove residual momen¬ 
tum from the largest remnant and in other cases these 
“jets” are symmetric and account for zero net momen¬ 
tum in the collision outcome. In order to determine the 
most accurate placement for the resolved fragments we 
analysed the debris distributions of the several hundred 
rubble-pile simulations used to create the EDACM colli- 
sion model presented in Leinhardt & Stewart (2012) and 



Figure 3. Post-collisional debris viewed face-on looking down on 
the collision in the u direction (see Fig. [4j|. Frames (a) and 
(b) show the result of a rubble-pile impact. Frames (c) and (d) 
show the result using the EDACM model. Both are integrated 
using PKDGRAV. The impact shown in (a) and (c) is a hit-and- 
run projectile disrupt (Vi = 50 m s 1 , oblique impact: b=0.9), 
between a km-sized target and a projectile a quarter of the mass. 
The impact shown in (b) and (d) involves the same projectile and 
target but is head-on and super-catastrophic (Vi =38 m s _1 . b = 
0 . 0 ). 


derived a series of rules for post-collision debris place¬ 
ment. Two very different collision debris fields from iso¬ 
lated rubble-pile collisions modelled with PKDGRAV are 
shown in Fig. [3jt-b. The EDACM generated version of 
these collisions is able to capture the general character 
of the debris field as shown in frames c-d. The details of 
the model implemented in EDACM are given below. 

Two angles are used to describe the distribution of 
fragment material 6 (0° < 9 < 360°), which deter¬ 
mines the direction of the debris swarms and cf> (—90° < 
<f> < 90°), which determines the degree of “jetiness” or 
“diskiness”. The fragment probability distribution (a bi- 
Gaussian distribution in 0 and a Gaussian distribution 
in <j>), 

-(e-ex) 2 -U -(e-p 2 ) 2 -A 

Pe<t> = A ei e e 2,7 i+Ag 2 e 2<T e e 2< ^, (8) 

is defined by two peaks in 6 located at 6 \, O 2 , and one at 
0 = 0° with thicknesses represented by the standard de¬ 
viations <jg and < 7 $ and relative amplitudes in 6 denoted 
as Ag 1 , Ag 2 . 

In order to determine the locations of all of the de¬ 
bris material let us first define a convenient collisional 
coordinate system u, v,w. The first basis vector, u is in 
the direction of the center joining vector, C = R p — R t , 
where R p and Rt are the positions of the projectile and 
target respectively (see Fig. [4]), thus, 

u=C/|C|. (9) 

The direction of the w basis vector is perpendicular to 
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Figure 4. Cartoon of the geometry of an impact between a pro¬ 
jectile (blue) and a target (red). The angle 0 indicates collisional 
debris axis. The value of <j> indicates whether the debris is dis¬ 
tributed in two jets or in a complete disk around the impact site. 


the impact velocity vector, Vi and C 

w = Vi/|V;| x u (10) 

and 

v = u x w. (11) 

Thus, the impact vector, V), is in the u — v plane. 

The debris axis is in the u — v plane and is defined by 
two peaks 


9i = 


70 b c . 


70 6 

t ^9° b 


• 6 ^ 6 crit 

90: 6 <c 6 crit 


e 2 = o 1 + 180, 


( 12 ) 

(13) 


with 9 = 0° in the u direction. If the collision was a non¬ 
grazing, head-on collision (6 = 0) the debris axis would 
be distributed about = 90° and 9 2 = 270°. Debris 
moves away from the impact site in two directions 180° 
apart with a narrow width that is relatively independent 
of collision parameters. Thus, the standard deviation of 
the debris distribution in 9 is set as a constant, 


ere = 10°. (14) 

The relative amplitudes of the two peaks in 9 depend 
strongly on impact parameter b and projectile to target 
mass ratio, /z, such that, 

~p~ = 1 — 6(1 — n). (15) 

Ag 2 

Thus, in all equal-mass collisions the debris is symmetric 
and independent of 6, however, for unequal mass impacts 
the debris can carry momentum away from the impact 
site. 

The peak of </> distribution, <fi = 0°, is in the v direc¬ 
tion. The width in <fi depends most prominently on the 
fraction of impacting mass, 

(70 = 145 a°, (16) 

where a is the projectile mass fraction involved in the 
collision (as defined in Eq. [2j . Note that the cf> distri¬ 
bution can be broad enough that the two debris “jets” 
actually touch creating a full disk (see Fig. §0- 


3. RESULTS 

The general evolution of the planetesimals is indepen¬ 
dent of the collision model (Fig. [5]) because most colli¬ 
sions in the relatively calm evolution scenario used in this 
work are accretion dominated and in addition, the veloc¬ 
ity dispersion is dominated by the gravity of the largest 
planetesimal bodies. This results in the largest planetes¬ 
imals following a similar growth history. The similarity 
seen here between the three collision models may not 
be seen in systems that undergo an energetic dynami¬ 
cal shake-up and/or have a large pert urber (for example, 
the Grand Tack scenario, Walsh et al. 2011). In any 
case, in the perturber-free scenario modelled here subtle 
differences amongst the different models are detectable 
in the evolution of the background planetesimal popula¬ 
tion which suffers more energetic collisions that result in 
fragmentation. As the largest embryos grow the RMS ve¬ 
locity of the background increases. The EDACM model 
allows fragmentation and provides both directional and 
velocity information about the collision remnant s, th us 
we can make a prediction of observable dust (see 4.1). 


3.1. Planetesimal Evolution 

Figure[5]shows snapshots of planetesimal evolution and 
embryo growth for three collision models: perfect merg¬ 
ing, RUBBLE and EDACM. The left hand panels show 
mass versus semi-major axis and right show eccentric¬ 
ity versus semi-major axis. The larger data points in 
these frames are used to indicate planetesi mals t hat have 
reached isolation mass (~ O.IM 0 ss 100m o , |Leinhardt & 
Richardson|2005 ). In all simulations the embryo popula- 
tion has begun to separate from the background popula¬ 
tion at small semi-major axis by 100,000 yrs (simulation 
time). By the end of the simulations all protoplanets 
are approaching a similar mass in dicating a transitio n 
from runaway to o l igarchic gro wth (Kokubo & Ida 2002 
Leinhardt & Richardson] 2005). 


As would be naively expected the perfect merging sim¬ 
ulation evolves the most quickly because it contains no 
fragmentation, then RUBBLE, and finally the EDACM 
simulation. In addition, there is no significant change to 
the growth modes due to different collision models. The 
only noticeable difference between the models in Fig. [5] is 
a slight increase in the evolution timescale due to less ef¬ 
ficient growth when fragmentation is included. This may 
get more pronounced with evolution time because the av¬ 
erage impact speeds increase and more of the collisions 
are erosive. None the less, once the difference in evolu¬ 
tion timescale is taken into account the three collision 
models qualitatively agree in eccentricity and inclination 
distributions, number, and mass of protoplanets. 

The RUBBLE and EDACM collision models allow the 
composition of all particles to be tracked. The compo¬ 
sition of the protoplanets (shown in Fig. [6] as pie plots) 
show some localised mixing. The pie plots are based on 
the composition histograms of each protoplanet which 
track where the mass accreted by each particle comes 
from. In general, most of the mass in each protoplanet 
comes from ±0.1 AU of the protoplanet’s current semi¬ 
major axis. This can also be seen by examining the 
colours of the protoplanets in Fig. [5] which indicate a 
mass-weighted average of the protoplanets composition 
(a mean of the pie plot). There is no significant differ- 
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Figure 5. Evolution of mass and eccentricity of resolved particles (planetesimals and protoplanets) versus semi-major axis for three 
collision models: merging, RUBBLE, and ED ACM. Protoplanets, objects that have reached isolation mass (100 times initial mass), are 
indicated with black filled circles (left) and coloured circles with error bars (right). Planet esimals, objects that are directly resolved but 
have not reached isolation mass, are shown as small open or coloured circles. Time increases from top to bottom. The simulation time 
for each row is indicated in black on far right. The red italicised time is the estimated effective time for un-inflated planetesimals with 
expansion parameter of one. Right - error bars show gravitational influence of protoplanets and extend 5 Hill radii on either side. The size 
of the protoplanets is scaled by radius. Rubble and ED ACM collision model columns also have colour derived from composition histograms 
associated with each particle. 
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Figure 6. Composition pie plots of each protoplanet from the RUBBLE and EDACM runs. The protoplanets are shown at their 
instantaneous semi-major axis. The size of the protoplanet is proportional to the object’s radius. Each wedge of a pie indicates the fraction 
of mass from the designated radial region shown by the colour bar at the bottom of the figure. 


ence in the degree of mixing between the two collision 
models. 

In all simulations, no matter what collision model is 
used, we find that the background planetesimal popula¬ 
tion is efficiently cleared out by the embryos. There is 
no significant period of time over the 1AU radial annu¬ 
lus in any simulation in which half the mass is in small 
planetesimals and half the mass is in large protoplanets 
(Fig-0) . Moreover, the difference in evolution timescale 
in the inner region of the annulus and the outer region 
of the annulus is noticeable even within 1 AU. The inner 
region of the annulus (inward of 0.8 AU) has been all 
but entirely cleared of planetesimals by the end of the 


simulations where as a considerable fraction of the mass 
if not the majority is in the planetesimal population in 
the outer regions of the annulus. These results bring 
into question argument s made in previous work such as, 


Goldreich et al. (2004), that a significant mass can be 
stored in the planetesimal population and thus reduce 
the orbital eccentricities of the oligarchies via dynamical 
friction at late times. 


3.2. Collisions 

Although most of the collisions in RUBBLE and 
EDACM simulations are accretion dominated very few 
of the collisions in the EDACM collision model actually 
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result in perfect merging (see black sections in each bar 
of Fig. [8]), instead the collision outcomes are divided be¬ 
tween partial accretion events (dark blue) and some type 
of hit-and-run (light blue - projectile intact, light green 
- projectile disrupted, and dark green - projectile catas¬ 
trophically disrupted). The fraction of partial accretion 
events remain relatively constant just above 40% while 
the fraction of hit-and-run events with an intact projec¬ 
tile (light blue) begins to decrease with time. This is 
due to the increase in the average impact speed as the 
largest bodies grow and gravitationally stir the rest of 
the objects. With an increase in speed even impacts 
with too much angular momentum to result in an accre¬ 
tion or merging event deliver enough energy to disrupt 
or partially disrupt the projectile. In addition to the de¬ 
crease in “perfect” hit-and-run events with time erosive 
events also begin to increase in number (yellow and red). 
Again this is due to the increased average speed of colli¬ 
sions due to gravitational stirring from growing planetary 
embryos. However, it is clear that the vast majority of 
all collisions are non-erosive (for the target) during all 
phases of evolution. 

The number of collisions decreases exponentially with 
time (white outlined histogram) as runaway growth 
quickly reduces the number of particles available for col¬ 
lisions. As a result, the collision type histogram be¬ 
comes more noisy with time as the low number of col¬ 
lisions reduces the statistical accuracy of each late time 
bin. The collisional breakdown shown in Fig. [8] shows 
slightly more hit-and-run projectile intact and slightly 
less catastrop hic disruption events than found in |Bon-| 
sor et al. (2015) which investigates the formation ofThe 
Earth using similar numerical simulations. The differ¬ 
ence in these results is due to a difference in initial size 
distribution and resolution limit and not to any subtle 
difference in the collision model. 

Figure [5] shows that varying the collision model does 
not significantly effect the growth phases of runaway or 
oligarchic growth due to the dominance of bouncing and 
accretion dominant collisions (Fig.pl). However, includ¬ 
ing a realistic fragmentation model does allow us to make 
observational predictions which would not otherwise be 
possible. 

4. DISCUSSION 

Although the direct numerical simulations presented 
here have a fixed resolution limit of ~ 500 km we can use 
the orbital information from the resolved planetesimals 
along with the mass in unresolved material (Fig. [9]) to 
help us make predictions of the location and intensify of 
dust visible during the formation of planetary embryos. 
In order to do this we follow the method described in 


Dobinson et al. (2013) which is briefly summarised in 
the section below. 


4.1. Synthetic Dust Images 

The first step in creating a synthetic image of dust is to 
determine the dust surface density. Predicting the dust 
surface density beyond the resolution limit of a numer¬ 
ical simulation has also been used in work invest igating 
the formation and detectability of debris disks (Booth 


et al. 


2009). However, in that case the authors assumed 


that the resolved planetesimals would never grow and 
provided only source material for dust. Thus, in the case 



0.5 0.75 1 1.25 1 . 50.5 0.75 1 1.25 1 . 50.5 0.75 1 1.25 1.5 

Semi-major Axis [AU] 


Figure 7. Fractional mass in protoplanets as a function of radial 
distance for the three collision models. The shaded histograms 
show the mass in protoplanets over a radial bin of 0.1 au normalised 
by the total mass in that bin. Time evolution is indicated by the 
different shades of grey increasing with time from dark to light. 
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Figure 8. Coloured histogram shows the percentage of planetesi- 
mal collision type per time bin for the extent of the simulation (left 
y-axis). Each time bin has a width of 2 X 10 5 yr in simulation time. 
The white outlined histogram indicates the number of collisions as 
a function of time (right y-axis). 

of debris disk formation the planetesimal mass could be 
used as a direct tracer for dust, namely, the more massive 
the planetesimal, the more dust generated in a collisional 
cascade. In the planet formation environment presented 
in this work the planetesimals are not in any form of 
steady state, instead they are growing and evolving due 
to collisions with each other. Therefore, we would expect 
dust to trace dynamical activity more directly than mass. 
In addition, we will also assume that all collisions pro¬ 
duce some dust even if the collision between the parent 
planetesimals results in a perfect merging event. The to¬ 
tal dust mass is determined by the mass flux of resolved 
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Figure 9. Evolution of unresolved debris in ED ACM simulation. 
The total mass of unresolved debris in units of initial planetesimal 
mass m 0 is show with open circles. Unresolved debris outside of the 
integration annulus, which lies between 0.5 and 1.5 AU, is shown 
with black crosses and is not available to be accreted by resolved 
planet esimals. 


which gives us a solid angle of 1.06 x 10~ 14 sr. This makes 
the flux per beam or resolution element 5.3 x 10~ 9 Jy. 
Observing using the full ALMA array at 353 GHz (850 
/rm) gives a sensitivity equal to the flux in 5.25x 10 s days. 
Obviously not practical, however, if the flux scaled with 
mass and the disk were 1000 times more massive such a 
disk would be detectable with a 12 hour observation. A 
disk with 1000 times more mass may not be that extreme. 
The protoplanetary disk shown in this work is only 3x3 
AU 2 which is over an order of magnitude smaller than the 
current extent of our own Solar System and almost two 
orders of magnitude smaller than the young protoplan¬ 
etary disk imaged by ALMA around HL Tau. The disk 
investigated in this work assumed a minimum mass solar 
nebula which is most likely an underestimate for our own 
solar system. The diversity of extrasolar planets and pro¬ 
toplanetary disks certainly suggest that accretion disks 
could be much more massive. In addition, the scenario 
presented in this paper involved no external perturbers 
such as a nearby giant planet or binary companion, which 
may significantly enhance planetesimal collision veloci¬ 
ties, dust production, and increase the possibility of ob¬ 
servational detection (see Dobinson et al. sub, Carter et 
al. in prep). Thus, if we observed the combined flux from 
a dust disk that contained 1000 times more dust mass it 
would be marginally detectable with ALM^J 


planetesimals on crossing orbits (see r eso lved planetesi¬ 
mal surface density in first row of Fig. 10). The mass in 
dust generated by predicted collisions with each planetes¬ 
imal is smeared over the orbit of the target planetesimal 
with mass distributed in equal time chunks, thus, highly 
eccentric orbits are not uniformly bright. The total mass 
in dust over the entire region is then scaled to the average 
total mass in unresolved debris close in time to each time 
step plotted (see Fig. 10 second row). The instantaneous 
dust mass is averageato reduce the noise in the dust 
surface density due to impact events at or near a given 
time step, which can be seen in the spikes and valleys in 
Fig. [9} 

A telescope would only observe a fraction of this total 
dust mass - namely, dust of a size similar to the observed 
wavelength. In addition, the brightness of the dust will 
depend strongly on the proximity of the dust to the cen¬ 
tral star, which is assumed t o be solar. Taking this into 
account the third row of Fig. [To] shows synthetic images 
produced using the radiative transfer code RADMC3D 
(Dullemond|2012|) at an ALMA observing wavelength of 
850 /rrn, assuming a maximum of one hundredth of the 
total unresolved debris mass at each time step ends up in 
dust between 0.1 and 1000 fj,m (roughly estimated from 
assuming a power-law size distribution with a slope of 
—3.5), where the lower limit on the dust size is the blow¬ 
out radius and the upper limit on the dust size is de¬ 
termined by the largest dust particle that significantly 
contributes to the sub-mm flux. The resulting images 
show a disk with a maximum flux of about 1 x 10 6 Jy/sr 
that both dims and becomes more ring-like as planetesi¬ 
mal evolution and protoplanet formation occur in ernest 
in the inner terrestrial region. 

What would be required to resolve such structures 
with ALMA? If we assume that the average brightness of 
5 x 10 5 Jy/sr and that we need to resolve disk structure 
on the scale of 0.2 AU at 10 pc then our beamsize is 0” .02 


5. CONCLUSION 

Large swaths of planet formation are currently unob¬ 
servable in a direct sense due to lack of emission from 
the evolving planetesimals. However, indirect observa¬ 
tion of planetesimal growth may be possible in the near 
future by observing collisionally generated dust with cut¬ 
ting edge observatories. In this paper we present our 
predictions for these observations using the efficient N- 
body code PKDGRAV to model runaway and oligarchic 
growth coupled with ED ACM, an empirically derived 
collision model that allows multiple collision outcomes 
including perfect merging, accretion, and erosion. We 
find that planetesimal growth and evolution could be in¬ 
directly observable with a fully functioning ALMA-like 
telescope via collisionally generated dust even in a calm 
system if the total dust production was large enough. 
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